# ANES models (Appendix)
# Load packages and install them if necessary (automatic)
need = c(
  "dplyr",           #data structuring/cleaning
  "estimatr",        #run lm_robust models
  "extrafont",       #plotting (fonts)
  "ggplot2",         #plotting
  "here",            #set wd
  "janitor",         #data structuring/cleaning
  "sandwich",        #robust se's
  "scales",          #convenience function to recode to 0-1
  "stringr",         #data structuring/cleaning
  "tidyr"            #data structuring/cleaning
  )     
##
have = need %in% rownames(installed.packages())
if ( any(!have) ) { install.packages( need[!have] ) }
pack <- lapply(need, library, character.only = TRUE)
rm(have,need,pack)

# Load fonts
extrafont::loadfonts()

# Start with setting a seed for replication purposes (just in case, also set in function)
set.seed(28)

# Optional but handy (turn off scientific notations)
options(scipen = 999) 

# Load functions for bootstraps and computing output
source("D:/dv_rep/Source/boot_functions.R")
# Load data
source("D:/dv_rep/Source/dat_anes.R")

#### Bootstrapped models ####

# Set N bootstraps to 1000
n_bootstraps <- 1000

# Bootstraps without controls (note: together, these bootstraps will take a considerable amount of time, up to several hours!)
boot_anesPolviews  <- perform_bootstrap(anes_Polviews, n_bootstraps, include_controls = FALSE)
boot_anesSwelfare  <- perform_bootstrap(anes_Swelfare, n_bootstraps, include_controls = FALSE)
boot_anesRacgender <- perform_bootstrap(anes_Racgender, n_bootstraps, include_controls = FALSE)
boot_anesCultmoral <- perform_bootstrap(anes_Cultmoral, n_bootstraps, include_controls = FALSE)
boot_anesBlkaid    <- perform_bootstrap(anes_Blkaid, n_bootstraps, include_controls = FALSE)

# Bootstraps with controls
cboot_anesPolviews  <- perform_bootstrap(anes_Polviews, n_bootstraps, include_controls = TRUE)
cboot_anesSwelfare  <- perform_bootstrap(anes_Swelfare, n_bootstraps, include_controls = TRUE)
cboot_anesRacgender <- perform_bootstrap(anes_Racgender, n_bootstraps, include_controls = TRUE)
cboot_anesCultmoral <- perform_bootstrap(anes_Cultmoral, n_bootstraps, include_controls = TRUE)
cboot_anesBlkaid    <- perform_bootstrap(anes_Blkaid, n_bootstraps, include_controls = TRUE)

# Not run: Run robustlm models with clustered SEs by year (these are not reported in the manuscript, but can be viewed for reference)
#source("D:/Source/robustlm_anes.R")

# Extract parameters models without controls
param_anesPolviews  <- extract(
  boot_anesPolviews, 
  data_frame  = anes_Polviews
)

param_anesSwelfare  <- extract(
  boot_anesSwelfare, 
  data_frame = anes_Swelfare
)

param_anesCmoral  <- extract(
  boot_anesCultmoral, 
  data_frame = anes_Cultmoral
)

param_anesRgender  <- extract(
  boot_anesRacgender, 
  data_frame = anes_Racgender
)

param_anesBlkaid  <- extract(
  boot_anesBlkaid, 
  data_frame = anes_Blkaid
)

# Models with control variables
cparam_anesPolviews  <- extract(
  cboot_anesPolviews, 
  data_frame = anes_Polviews
)

cparam_anesSwelfare  <- extract(
  cboot_anesSwelfare, 
  data_frame = anes_Swelfare
)

cparam_anesCmoral  <- extract(
  cboot_anesCultmoral, 
  data_frame = anes_Cultmoral
)

cparam_anesRgender  <- extract(
  cboot_anesRacgender, 
  data_frame = anes_Racgender
)

cparam_anesBlkaid  <- extract(
  cboot_anesBlkaid, 
  data_frame = anes_Blkaid
)

# Make data frames for each model
anes_BootPolviews     <- make_Tab(param_anesPolviews)
anes_BootPolviews$var <- "Ideological self-placement"
anes_BootSwelfare     <- make_Tab(param_anesSwelfare)
anes_BootSwelfare$var <- "Social welfare"
anes_BootRgender      <- make_Tab(param_anesRgender)
anes_BootRgender$var  <- "Race and gender"
anes_BootCmoral       <- make_Tab(param_anesCmoral)
anes_BootCmoral$var   <- "Culture and morality"
anes_BootBlkaid       <- make_Tab(param_anesBlkaid)
anes_BootBlkaid$var   <- "Aid to Blacks"

# Models with control variables
canes_BootPolviews     <- make_Tab(cparam_anesPolviews)
canes_BootPolviews$var <- "Ideological self-placement"
canes_BootSwelfare     <- make_Tab(cparam_anesSwelfare)
canes_BootSwelfare$var <- "Social welfare"
canes_BootRgender      <- make_Tab(cparam_anesRgender)
canes_BootRgender$var  <- "Race and gender"
canes_BootCmoral       <- make_Tab(cparam_anesCmoral)
canes_BootCmoral$var   <- "Culture and morality"
canes_BootBlkaid       <- make_Tab(cparam_anesBlkaid)
canes_BootBlkaid$var   <- "Aid to Blacks"

# Combine to a single data frame
anes_BootDat <- plyr::rbind.fill(
  anes_BootPolviews,
  anes_BootSwelfare,
  anes_BootRgender,
  anes_BootCmoral,
  anes_BootBlkaid
) 

# Add column to specify that these models are without control variables
anes_BootDat$controls <- paste("Without")

# Models with controls
canes_BootDat <- plyr::rbind.fill(
  canes_BootPolviews,
  canes_BootSwelfare,
  canes_BootRgender,
  canes_BootCmoral,
  canes_BootBlkaid
) 

# Repeat
canes_BootDat$controls <- paste("With")

# Combine to single data frame
anes_Boot <- plyr::rbind.fill(
  anes_BootDat,
  canes_BootDat
)

# Data structuring: add columns for indexing
anes_Boot <- anes_Boot %>%
  mutate(term = gsub(":_C", "", term),
         term = case_when(
           term == "_C" ~ "Variable",
           term == ":white" ~ ":White",
           term == ":college" ~ ":College", 
           term == ":protcattend" ~ ":Protestant",
           TRUE ~ term
         )) %>%
  mutate(group = case_when(
    startsWith(term, "Cohort") ~ "Cohort",
    startsWith(term, "Age") ~ "Age",
    startsWith(term, "Year") ~ "Year",
    startsWith(term, ":Cohort") ~ "Cohort",
    startsWith(term, ":Age") ~ "Age",
    startsWith(term, ":Year") ~ "Year",
    startsWith(term, ":_C:Cohort") ~ "Cohort",
    startsWith(term, ":_C:Age") ~ "Age",
    startsWith(term, ":_C:Year") ~ "Year",
    startsWith(term, "White") ~ "White",
    startsWith(term, "College") ~ "College",
    startsWith(term, "Year") ~ "Protestant",
    startsWith(term, ":White") ~ "White",
    startsWith(term, ":College") ~ "College",
    startsWith(term, ":Protestant") ~ "Protestant",
    TRUE ~ term
  )) %>%
  mutate(boot = case_when(
    startsWith(model, "Boot") ~ "Boot",
    TRUE ~ "OLS"
  )) %>%
  mutate(term = gsub("Cohort |Age |Year ", "", term),
         facet = gsub("Boot |Original ", "", model)) 

# Add proportion of respondents for each group to data frame (we'll need this to compute values for plotting)
# This proportion column is called Freq (frequency respondents)
prop_polviews  <- prop(anes_Polviews, "Ideological self-placement", gen = "genPew")
prop_swelfare  <- prop(anes_Swelfare, "Social welfare", gen = "genPew")
prop_rgender   <- prop(anes_Racgender, "Race and gender", gen = "genPew")
prop_cultmoral <- prop(anes_Cultmoral, "Culture and morality", gen = "genPew")
prop_blkaid    <- prop(anes_Blkaid, "Aid to Blacks", gen = "genPew")

prop_merge     <- plyr::rbind.fill(
  prop_polviews,
  prop_swelfare,
  prop_rgender,
  prop_cultmoral,
  prop_blkaid
)

anes_Boot <- anes_Boot %>% left_join(prop_merge, by = c("term", "var"))

# Save for future use
save(anes_Boot, file = "D:/dv_rep/Datasets/anes_Boot.Rdata")

# Additional data structuring for p-values
# Warning refers to the fact that we don't have a standard error for R2 for the OLS models, which is true
table_ANES <- anes_Boot %>% 
  mutate(
    p_star = case_when(
    p < 0.001 ~ "***",
    p < 0.01  ~ "**",
    p < 0.05  ~ "*",
    TRUE ~ ""
  )) %>%
  mutate(across(c(coef, se), ~ as.numeric(.))) %>%
  mutate(across(c(r2, r2_dev, adj_r2), ~ as.numeric(format_fixed_decimals(., decimals = 3)))) %>%
  mutate(coef_p = paste0(format_fixed_decimals(coef, 2), " (", format_fixed_decimals(se, 2), ")", p_star)
  )

# Convert to wide format
tab_Wide <- table_ANES %>%
  select(term, controls, coef_p, facet, boot, var) %>%
  pivot_wider(., names_from = var, values_from = coef_p) %>%
  janitor::clean_names() %>%
  select(
    term,
    facet,
    boot,
    controls,
    ideological_self_placement,
    social_welfare,
    race_and_gender, 
    aid_to_blacks,
    culture_and_morality
  ) %>%
  as.data.frame()

# Extract parameters of interest and add as new rows
rows_Tab <- table_ANES %>%
  select(controls, facet, boot, var, r2, r2_dev, adj_r2, N) %>%
  distinct %>%
  pivot_longer(cols = c(r2, r2_dev, adj_r2, N),
               names_to = "term", 
               values_to = "coef_p") %>%
  mutate(coef_p = ifelse(term %in% c("r2", "adj_2", "r2_dev"), format_fixed_decimals(coef_p, 3), coef_p)) %>%
  mutate(coef_p = as.character(coef_p)) %>%
  mutate(coef_p = case_when(term == "r2" & boot == "Boot" ~ paste0(coef_p, " (", lead(coef_p), ")"), TRUE ~ coef_p)) %>%
  filter(term != "r2_dev") %>%
  pivot_wider(., names_from = var, values_from = coef_p) %>%
  janitor::clean_names()
  
# Create a list of unique groups in tab_Wide
unique_groups <- unique(tab_Wide[c("controls", "facet", "boot")])

# Initialize an empty data frame to store the combined result
combined_df <- data.frame()

# Iterate over unique groups and append rows_Tab after each group in tab_Wide
for (i in 1:nrow(unique_groups)) {
  group <- unique_groups[i, ]
  group_rows_Tab <- rows_Tab %>%
    filter(controls == group$controls, facet == group$facet, boot == group$boot)
  
  group_tab_Wide <- tab_Wide %>%
    filter(controls == group$controls, facet == group$facet, boot == group$boot)
  
  combined_group <- bind_rows(group_tab_Wide, group_rows_Tab)
  combined_df <- bind_rows(combined_df, combined_group)
}
# Mean values
meanTab <- combined_df %>%
  filter(facet == "Mean" & boot == "Boot") %>%
  select(!boot & !facet)

meanTab <- meanTab %>% 
  pivot_wider(., names_from = controls, values_from = c(
    "ideological_self_placement",
    "social_welfare",
    "race_and_gender",
    "aid_to_blacks",
    "culture_and_morality"
  )) %>%
  mutate_all(~ replace_na(., "")) 
write.csv(meanTab, file = "meananes_Boot.csv", row.names = FALSE)
# Ideological divergence
divTab <- combined_df %>%
  filter(facet == "Divergence" & boot == "Boot") %>%
  select(!boot & !facet)
divTab <- divTab %>% 
  pivot_wider(., names_from = controls, values_from = c(
    "ideological_self_placement",
    "social_welfare",
    "race_and_gender",
    "aid_to_blacks",
    "culture_and_morality"
  )) %>%
  mutate_all(~ replace_na(., "")) 
write.csv(divTab, file = "divanes_Boot.csv", row.names = FALSE)
# Sorting
sortTab <- combined_df %>%
  filter(facet == "Sorting" & boot == "Boot") %>%
  select(!boot & !facet)
sortTab <- sortTab %>% 
  pivot_wider(., names_from = controls, values_from = c(
    "ideological_self_placement",
    "social_welfare",
    "race_and_gender",
    "aid_to_blacks",
    "culture_and_morality"
  )) %>%
  mutate_all(~ replace_na(., "")) 
# Write table to csv
write.csv(sortTab, file = "sortanes_Boot.csv", row.names = FALSE)
